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Molybdenum disulfide (M0S2) is a layered semiconductor which has become very important recently as 
an emerging electronic device material. Being an intrinsic semiconductor the two-dimensional M0S2 has 
major advantages as the channel material in field-effect transistors. In this work we determine the electronic 
structure of M0S2 with the highly accurate screened hybrid functional within the density functional theory 
(DFT) including the spin-orbit coupling. Using the DFT electronic structures as target, we have developed a 
single generic tight-binding (TB) model that accurately produces the electronic structures for three different 
forms of M0S2 - bulk, bilayer and monolayer. Our TB model is based on the Slater-Koster method with 
non-orthogonal sp^d^ orbitals, nearest-neighbor interactions and spin-orbit coupling. The TB model is useful 
for atomistic modeling of quantum transport in M0S2 based electronic devices. 
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Molybdenum disulfide (M0S2) belongs to a family 
of layered transition metal dichalcogenides (TMDC) in 
which the layers are held together by weak van der Waals 
forces, and it can be exfoliated mechanically to a sin- 
gle layer thickness. In its bulk form M0S2 is an indi- 
rect band gap semiconductor which turns into a direct 
band gap semiconductor for monolayer structure.^ This 
intrinsic semiconducting nature of M0S2 is a major ad- 
vantage over graphene (which has no intrinsic band gap) 
as a two-dimensional (2D) channel material in field-effect 
transistors (FET). Indeed, FET devices based on M0S2 
monolayer and bilayer have already been fabricated in 
the experimental labs and demonstrated to have use- 
ful device performances.^'^ Several theoretical studies of 
M0S2 FET devices^^ based on simplified description of 
the electronic structures within an effective mass model 
have also been reported recently. 

Due to the complex nature of the layered TMDC mate- 
rials and the great interests in the electronics community 
for its applications in emerging devices, a more reliable, 
accurate, and atomistic treatment of the electronic struc- 
tures of TMDC is desired. Since ah initio models are 
computationally expensive and often intractable for re- 
alistic device structures having large number of atoms, 
a widely applied intermediary option is the tight-binding 
(TB) model. TB models are atomistic full-band yet com- 
putationally efficient due to the small number of basis 
orbitals and the rigid nature of the Hamiltonian. Prop- 
erly parameterized TB models have been extensively ap- 
plied to predict reliable results for a diverse range of 
devices including the resonant tunneling diodes,^ quan- 
tum dots,^ Si nanowireFET,^ and compound semicon- 
ductor heterostructures.^^ However, for the very excit- 
ing 2D M0S2 and other TMDC materials, accurate TB 
models for electronic structures has not yet been devel- 



oped. It is the purpose of this work to fill this void, and 
the resulting TB model should be extremely useful for 
predicting quantum transport properties of M0S2 based 
devices. 

In particular, we have developed a generic TB model 
with non-orthogonal sp^d^ orbitals including the nearest- 
neighbour interactions and spin-orbit coupling, that ac- 
curately determines the electronic structures of bulk, 
monolayer and bilayer M0S2. A reliable and accurate 
target band structure is the primary requirement for a 
successful tight-binding modeling. In this study, for the 
calculations of the target band structures we employ den- 
sity functional theory (DFT) with the screened hybrid 
functional of Heyd, Scuseria, and Ernzerhof^^ (HSE) that 
has been shown in the literature^^^^ to produce accu- 
rate band gaps and reasonable effective masses for a wide 
range of semiconductors. Our calculated band structures 
of M0S2 from the self-consistent DFT-HSE shows excel- 
lent agreement with the available experimental data and 
compare well with other theoretical studies. ^^^^ Even 
though an accurate target band structure is obtained by 
DFT-HSE, it is nontrivial to develop a generic TB model 
which can accurately capture electronic structures and ef- 
fective masses of not only bulk M0S2 but also its double- 
layer and single-layer forms. For the layered TMDC 
structures there is a thickness dependent interlayer inter- 
action that significantly affects the electronic properties, 
the TB model must capture these effects so as to be ap- 
plicable to devices made of different layer thicknesses. As 
shown below, our generic TB model accurately produces 
the band structures for three different forms of M0S2: 
bulk, monolayer and bilayer. The reported TB model 
can be easily extended to other M0S2 structures of higher 
number of layers as well as other layered TMDC materi- 
als such as WS2, WSe2, and MoSe2. The accuracy of our 
TB model is vahdated by comparing the band gaps and 
effective masses with our calculated ah initio results. 
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FIG. 1. (color online) Band structures for monolayer (IL), 
bilayer (2L), and bulk M0S2 using screened hybrid density 
functional theory (blue lines) and our tight-binding fitting 
(red dots) . The zero in the energy axis is set at the Fermi level 
as shown by the dashed line. Spin-orbit coupling is included 
in the calculations self-consistently. 



In our study, all the DFT-HSE calculations are per- 
formed using the projector augmented- wave (PAW) pseu- 
dopotential plane-wave method as implemented in the 
VASP software package. ^^ A Monkhorst-Pack scheme is 
adopted for the sampling of the Brillouin zone with k- 
point meshes of 7 x 7 x 2 for the bulk and 9 x 9 x 1 for 
the two-dimensional (2D) structures. An energy cutoff 
of 280 eV is used in a plane wave basis set. We also in- 
cluded the spin-orbit coupling self-consistently. Due to 
the presence of van der Waals interaction which deter- 
mines the interlayer distance in M0S2, geometry opti- 
mization becomes quite tricky and complicated. With- 
out the van der Waals intereaction the lattice parameter 
c (vertical length of the unit cell) is overestimated by a 
large margin. However, it is not yet possible to include 
van der Waals interaction in the HSE calculations, and it 
is an area of active research. We therefore optimize the 
bulk crystal structure using the Perdew-Burke-Ernzerhof 
(PBE) functional^^ for which the van der Waals interac- 
tion can be included explicitly by adopting the DFT-D2 
method. ^^ We then use this optimized structure for the 
HSE calculations. The values of our optimized lattice pa- 
rameters for bulk M0S2 are: a = 3.179 A and c = 12.729 
A with a layer thickness of 3.135 A (S-S vertical distance). 
Although the value of c parameter improves significantly 
with the DFT-D2 method, it is still overestimated by 
around 3.6% compared to the experimental value. 



It is a well-known fact that DFT within the local den- 
sity approximation (LDA) or the generalized gradient ap- 
proximation (GGA) underestimates band gaps of semi- 
conductors due to factors such as the self interactions. 
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FIG. 2. (color online) Contour plots of the energies around 
the valence band maximum (VBM) and conduction band min- 
imum (CBM) for monolayer M0S2 from screened hybrid den- 
sity functional theory (blue solid lines) and our tight-binding 
model (red dashed lines). The contour lines are associated 
with A^ = ± ^ 
the Ah 



3Q , ^ ^^, zb^ eV, respectively and the Akx and 
are in the unit of — , where a is the unit cell length. 



This error can be minimized in the HSE scheme of hybrid 
functional method^^ by incorporating 25% short-range 
(SR) exact Hartree-Fock (HF) exchange with the PBE 
functional. A screening parameter /i defines the range 
separation, and is empirically set to 0.2 A~^ (HSE06 
scheme) for both the HF and PBE parts. In our calcu- 
lations for M0S2 we observe that the band gap energies 
vary significantly with the screening parameter for both 
the bulk and the 2D structures. We therefore decided to 
adjust the value of /i to fit the experimental band gap, 
and an optimize value for /i is found to be 0.4 A~^ which 
corresponds to an interaction range of 7r//i = 7.85 A for 
the SR nonlocal exchange. Note that a single value of /i 
is used for all the three structures of M0S2: bulk, double- 
layer and single-layer, namely the HSE results presented 
here were all obtained using /i = 0.4 A~^. 

Our calculated band structure results for M0S2 are 
plotted in Fig. 1. We notice that the monolayer M0S2 
has a direct band gap at the K point. However, with the 
addition of just one extra layer (i.e. for a bilayer struc- 
ture) an indirect band gap opens up due to the presence 
of interlayer interactions. Although both bulk and bi- 
layer M0S2 have indirect band gaps the conduction band 
minima (CBM) are located at different points: for bulk it 
is at S point (midpoint of F and K points) and for bilayer 
it is at K point. Similar behavior is observed in a previous 
theoretical study. ^^ For monolayer M0S2 the splitting of 
the valence band maximum (VBM) at K point is solely 
due to the spin-orbit coupling (SOC) whereas a combina- 
tion of SOC and interlayer interactions is responsible for 
the VBM splitting in bulk and bilayer M0S2. The split- 
ting of CBM due to SOC is minimal (around 5 meV). The 
calculated band gap energies for different transitions are 
presented in Table I. Our DFT-HSE results show good 
agreement with the experimental data. The values of ef- 
fective masses at different band edges along different di- 
rections as presented in Table II also compare well with 
other available theoretical studies. ^^'^^ 

With the calculated DFT-HSE band structures as tar- 
get, we employed the Nanoskif^^ software package to ob- 



TABLE I. Band gap energies obtained by DFT-HSE and our 
TB model. The fifth column is the deviation between the 
HSE and the TB values. Experimental data is shown in the 
sixth column, taken from Ref. 21 (for bulk) and Ref. 1 (for 
monolayer and bilayer). All the energies are in the unit of 
eV. Subscripts v and c stand for valence band and conduc- 
tion band, respectively. The splitting of the valence band 
maximum at K point is given by K^i (top) and Ky2 (bot- 
tom), whereas S is the midpoint of the line joining the F and 
the K points. 



Structure 


Gap 


HSE 


TB 


Deviation 


Exp. 






(target) 


(fitted) 


(%) 




Monolayer 


K^i to Kc 


1.786 


1.805 


1.06 


1.90 




K^2 to Kc 


1.974 


1.969 


0.24 


2.05 


Bilayer 


r^ to Kc 


1.480 


1.516 


2.41 


1.60 




K^i to Kc 


1.779 


1.792 


0.76 


1.88 




K^2 to Kc 


1.980 


1.987 


0.35 


2.05 


Bulk 


r^ to Sc 


1.328 


1.331 


0.22 


1.29 




K^i to Kc 


1.776 


1.749 


1.46 


1.88 




K^2 to Kc 


1.960 


2.009 


2.46 


2.06 



tain a set of TB parameters for the on-site energies, the 
Slater-Koster energy integrals, overlap integrals, ^^ and 
the spin-orbit splitting. ^^ In our TB model we adopted 
a non-orthogonal basis set of sp^d^ orbitals considering 
only the nearest neighbour interactions. The effects of 
spin-orbit coupling is included through a split-off energy 



term. A set of 96 parameters are optimized through the 
automated process built in Nanoskif.^^ The root-mean- 
square deviation of the fitting is within 25 meV. The op- 
timized TB parameters are listed in Table IIL Note that 
this single set of TB parameters is capable of producing 
accurate band structures of all three different forms of 
M0S2. 

The band structures obtained from our TB model are 
plotted in Fig. 1 along with those obtained from DFT- 
HSE. The agreement in energies around both VBM and 
CBM for all three structures of M0S2 is excellent. Fig. 2 
shows the energy contours around the valence band maxi- 
mum (VBM) and the conduction band minimum (CBM) 
for the monolayer M0S2. Again, we observe excellent 
agreement in the shape of the energy contours from DFT- 
HSE and TB fitting. The band gap energies for different 
transitions and the effective mass values at different sym- 
metry points obtained from the TB model are presented 
in Table I and Table H respectively. For the band gap en- 
ergies a fitting accuracy of less than 2.5% is achieved. On 
the other hand, for the effective masses, in most cases the 
deviations between the DFT-HSE and the TB values are 
less than 10% which is quite acceptable. In four cases 
(excluding the monolayer structure) the deviations are 
higher with the worst situation at ~ 27.9%. Given that 
just a single set of TB parameters is used to produce 
the band structurues of all three different forms of the 
material, this level of quantitative consistency is rather 
satisfactory. 



TABLE II. Values of effective masses at various band edges in the unit of free electron mass (mo) calculated using the HSE 
method and our tight-binding model. The subscripts / and t refer to the masses calculated at the point along the longitudinal 
and the transversal directions of the line connecting the F point and that point, respectively. 
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Point 


HSE 


TB 


Deviation 
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TB 


Deviation 






(target) 


(fitted) 


(%) 


(target) 


(fitted) 


(%) 


Monolayer 


K^ 


0.407 


0.430 


4.58 


0.485 


0.463 


4.54 




Kt 


0.404 


0.426 


5.45 


0.480 


0.458 


4.58 


Bilayer 


F 








1.039 


1.321 


27.91 




Ki 


0.430 


0.457 


6.28 










Kt 


0.423 


0.454 


7.33 








Bulk 


F 








0.785 


0.917 


16.82 




^i 


0.574 


0.712 


24.04 










^t 


0.819 


0.999 


21.98 









In conclusion, we have developed a generic TB model 
for accurately calculating the band structures of bulk, 
monolayer and bilayer M0S2. Our TB model is based on 
the Slater-Koster method. For the optimization of the 
TB parameters, accurate target band structures are ob- 
tained using the screened hybrid DFT method. The ac- 



I 

curacy of our TB model is verified by comparing the band 
gaps for different transitions and the effective masses at 
different band edges against the ab initio band structures. 
One main feature of our TB model is that with only one 
set of parameters it can reproduce the band structures 
of M0S2 of different structural configurations: from bulk 



TABLE III. Tight-binding parameters for M0S2 using non-orthogonal model with sp'^d^ orbitals, nearest-neighbour interactions, 
and spin-orbit coupling: on-site energies (E), spin-orbit splitting (A), Slater-Koster energy integrals {Ei for intra-layer and E2 
for inter- layer interaction) and overlap integrals (Oi for intra-layer and O2 for inter-layer interaction). The energies are in the 
unit of eV. 





Es 
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Ed 


Aso 










s 


7.6595 


-2.1537 


8.7689 


0.2129 










Mo 


5.5994 


6.7128 


2.6429 


1.0675 












Ei{^Mo) 


^i(S,S) ^2(S,S) E, 


(Mo,Mo) 


Oi(S,Mo) 


Oi(S,S) 


02(S,S) 0, 


.(Mo,Mo) 


ssa 


-0.0917 


0.3093 


0.3207 


0.1768 


0.0294 


-0.0532 


-0.1430 


-0.0575 


spa 


0.6656 


-0.9210 


-0.1302 


1.0910 


0.1042 


0.0240 


0.0196 


0.0057 


psa 


-1.6515 








0.1765 








ppa 


1.4008 


0.7132 


0.7053 


-0.3842 


-0.1865 


0.0478 


-0.0486 


0.0296 


PP7T 


-0.4812 


-0.1920 


-0.0980 


0.5203 


0.0303 


-0.0104 


0.0117 


0.0946 


sda 


0.2177 


-0.2016 


0.1164 


-0.5635 


-0.0480 


0.0946 


0.0297 


-0.1082 


dsa 


-1.0654 








-0.1432 








pda 


-2.8732 


-0.5204 


-0.0334 


-0.2316 


0.0942 


0.0724 


-0.0087 


0.0212 


dpa 


2.1898 








0.2002 








pdn 


0.7739 


-0.1203 


-0.0370 


0.0582 


0.0132 


0.0772 


-0.0031 


-0.0448 


dpTT 


-1.9408 








-0.2435 








dda 


-3.1425 


0.8347 


-0.2300 


0.3602 


0.0273 


0.1849 


0.0060 


-0.0216 


ddn 


2.4975 


0.7434 


0.0050 


0.0432 


0.1940 


-0.0429 


-0.0378 


-0.0285 


dd6 


-0.3703 


-0.1919 


-0.1104 


0.1008 


0.1261 


-0.0333 


0.0007 


0.0432 



to 2D structures. Our tight-binding model can be easily 
extended to other TMDC materials that show electronic 
characteristics similar to M0S2. The TB model reported 
in this paper is useful for simulations of quantum trans- 
port in nanoelectronic devices based on the TMDC ma- 
terials. 
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